home *** CD-ROM | disk | FTP | other *** search
/ NeXT Education Software Sampler 1992 Fall / NeXT Education Software Sampler 1992 Fall.iso / Programming / Source / 2DLab / vregion.c < prev   
Encoding:
C/C++ Source or Header  |  1992-01-14  |  14.1 KB  |  570 lines
  1. #include "headers.h"
  2.  
  3. #define LBND 0
  4. #define BBND 1
  5. #define RBND 2
  6. #define TBND 3
  7.  
  8. /* left, right, top bottom */
  9. #define NBOUND(x,y) (x >= Pxmax ? RBND : (x <= Pxmin ? LBND : (y >= Pymax ? TBND : (y <= Pymin ? BBND : -1))))
  10.  
  11. #define BBOX_THERE() ((pxmin == txmin) && (pxmax == txmax) && (pymin == tymin) && (pymax && pymax))
  12. #define BBOX_NOT_THERE() (!BBOX_THERE())
  13.  
  14. static float Pxmin, Pxmax, Pymin, Pymax;    /* floating point tolerance values */
  15. static float pxmin = 1E10, pxmax = -1E10, pymin = 1E10, pymax = -1E10;
  16. /* space for theta after x*y* */
  17. static float xminymin[3], xminymax[3], xmaxymax[3], xmaxymin[3];
  18. /* these point to the x*y* */
  19. static VERTTHETA *corner[4];
  20.  
  21. static POINT v, v1, v2, v3, v4;
  22.  
  23. static int numverts = 0, numvedges = 0, numtris;
  24.  
  25. /* sites[i].coord.x,y defined and sorted in main() */
  26.  
  27. VERT   GLsites[MAXVERTS];
  28. static VERT verts[MAXVERTS];
  29. static EDGE vedges[MAXEDGES];
  30. static TRI  tris[MAXTRIS];
  31. static int  inverse[MAXVERTS];
  32.  
  33. float randr(low,high)
  34. float low, high;
  35. {
  36.     float value;
  37.  
  38.     value = low + (rand()/((1<<15)-1.0))*(high-low);
  39.     return( value );
  40. }
  41.  
  42. out_site(s)
  43. struct Site *s;
  44. {
  45. }
  46.  
  47. out_bisector(e)
  48. struct Edge *e;
  49. {
  50. }
  51.  
  52. out_ep(e)
  53. struct Edge *e;
  54. {
  55.     if(!triangulate)
  56.     clip_line(e);
  57.  
  58.     if(debug)    {
  59.     printf("out_ep(): edge %d", e->edgenbr);
  60.     printf(" ep %d", e->ep[le] != (struct Site *)NULL ? e->ep[le]->sitenbr : -1);
  61.     printf(" ep %d", e->ep[re] != (struct Site *)NULL ? e->ep[re]->sitenbr : -1);
  62.     printf(" reg %d", e->reg[le] != (struct Site *)NULL ? e->reg[le]->sitenbr : -1);
  63.     printf(" reg %d\n", e->reg[re] != (struct Site *)NULL ? e->reg[re]->sitenbr : -1);
  64.     }
  65. }
  66.  
  67. out_vertex(v)
  68. struct Site *v;
  69. {
  70.     if (!triangulate)    {
  71.     verts[numverts].x = v->coord.x;
  72.     verts[numverts].y = v->coord.y;
  73.     if (numverts < MAXVERTS)
  74.         numverts++;
  75.     else    {
  76.         fprintf (stderr, "\nvert list overflow!");
  77.         exit (-1);
  78.         }
  79.     }
  80.  
  81.     if(debug)
  82.     printf("vertex(%d) at %10.7f %10.7f\n", v->sitenbr, v->coord.x, v->coord.y);
  83. }
  84.  
  85. out_triple(s1, s2, s3)
  86. struct Site *s1, *s2, *s3;
  87. {
  88.     if(triangulate) {
  89.     tris[numtris].v1 = (VERT *)&(sites[s1->sitenbr].coord);
  90.     tris[numtris].v2 = (VERT *)&(sites[s2->sitenbr].coord);
  91.     tris[numtris].v3 = (VERT *)&(sites[s3->sitenbr].coord);
  92.  
  93.     if (numtris < MAXTRIS)
  94.         numtris++;
  95.     else    {
  96.         fprintf (stderr, "out_triple(): triangle list overflow!\n");
  97.         exit (-1);
  98.         }
  99.     }
  100. }
  101.  
  102. vdinit()
  103. {
  104.     numverts = numvedges = numtris = 0;
  105. }
  106.  
  107. /* here, want to copy the gl sites INTO the voronoi array */
  108. /* gl sites are written exactly once at beginning of time */
  109. bboxinit()
  110. {
  111. int k;
  112. float dx, dy, x, y;
  113.  
  114.     /* get tight bounding box */
  115.     xmin =  1e10;
  116.     xmax = -1e10;
  117.  
  118.     ymin =  1e10;
  119.     ymax = -1e10;
  120.  
  121.     for(k = 0; k < nsites; k++) {    
  122.     x = sites[k].coord.x = GLsites[k].x;
  123.     y = sites[k].coord.y = GLsites[k].y;
  124.  
  125.         sites[k].refcnt = 0;
  126.     sites[k].sitenbr = k;
  127.  
  128.     if (x < xmin) xmin = x;
  129.     if (y < ymin) ymin = y;
  130.  
  131.     if (x > xmax) xmax = x;
  132.     if (y > ymax) ymax = y;
  133.     }
  134.  
  135.     /* NOW: xmin, ymin, xmax, ymax EXACT, as required by voronoi() */
  136.     /* we shall fool with pxmin, pymin, pxmax, pymax, as used by clip() */
  137.  
  138. #define EPSILON 1.0
  139.  
  140.     /* compute 'loose' bounding box */
  141.     dx = xmax - xmin;
  142.     dx = MAX (dx, EPSILON);
  143.  
  144.     dy = ymax - ymin;
  145.     dy = MAX (dy, EPSILON);
  146.  
  147.     pxmin = xmin - dx * 0.25;
  148.     pymin = ymin - dy * 0.25;
  149.  
  150.     pxmax = xmax + dx * 0.25;
  151.     pymax = ymax + dy * 0.25;
  152.  
  153. #ifdef STUPID_C_COMPILER
  154.     printf ("/* xmin, ymin %10.7f %10.7f; xmax, ymax %10.7f %10.7f; */\n", xmin, ymin, xmax, ymax);
  155.     printf ("/* pxmin, pymin %10.7f %10.7f; pxmax, pymax %10.7f %10.7f; crad %10.7f; */\n", pxmin, pymin, pxmax, pymax, cradius);
  156. #endif
  157. }
  158.  
  159. int clip_line(e)
  160. struct Edge *e;
  161. {
  162. struct Site *s1, *s2;
  163. struct Site *r1, *r2;
  164. float x1,x2,y1,y2;
  165.  
  166.     if(e->a == 1.0 && e->b >= 0.0) {    
  167.     s1 = e->ep[1];
  168.     s2 = e->ep[0];
  169.  
  170.     r1 = e->reg[1];
  171.     r2 = e->reg[0];
  172.     } else {    
  173.     s1 = e->ep[0];
  174.     s2 = e->ep[1];
  175.  
  176.     r1 = e->reg[0];
  177.     r2 = e->reg[1];
  178.     }
  179.  
  180.     if(e->a == 1.0) {
  181.     y1 = pymin;
  182.     if (s1!=(struct Site *)NULL && s1->coord.y > pymin)
  183.         y1 = s1->coord.y;
  184.     if(y1>pymax)
  185.         return;
  186.     x1 = e->c - e->b * y1;
  187.     y2 = pymax;
  188.     if (s2!=(struct Site *)NULL && s2->coord.y < pymax)
  189.         y2 = s2->coord.y;
  190.     if(y2<pymin)
  191.         return(0);
  192.     x2 = e->c - e->b * y2;
  193.     if ((x1> pxmax && x2>pxmax) || (x1<pxmin && x2<pxmin)) 
  194.         return;
  195.     if(x1> pxmax) {
  196.         x1 = pxmax;
  197.         y1 = (e->c - x1)/e->b;
  198.     }
  199.     if(x1<pxmin) {
  200.         x1 = pxmin;
  201.         y1 = (e->c - x1)/e->b;
  202.     }
  203.     if(x2>pxmax) {
  204.         x2 = pxmax;
  205.         y2 = (e->c - x2)/e->b;
  206.     }
  207.     if(x2<pxmin) {
  208.         x2 = pxmin;
  209.         y2 = (e->c - x2)/e->b;
  210.     }
  211.     } else {
  212.     x1 = pxmin;
  213.     if (s1!=(struct Site *)NULL && s1->coord.x > pxmin) 
  214.         x1 = s1->coord.x;
  215.     if(x1>pxmax)
  216.         return(0);
  217.     y1 = e->c - e->a * x1;
  218.     x2 = pxmax;
  219.     if (s2!=(struct Site *)NULL && s2->coord.x < pxmax)
  220.         x2 = s2->coord.x;
  221.     if(x2<pxmin)
  222.         return(0);
  223.     y2 = e->c - e->a * x2;
  224.     if ((y1> pymax && y2>pymax) || (y1<pymin && y2<pymin))
  225.         return(0);
  226.     if(y1> pymax) {
  227.         y1 = pymax;
  228.         x1 = (e->c - y1)/e->a;
  229.     }
  230.     if(y1<pymin) {
  231.         y1 = pymin;
  232.         x1 = (e->c - y1)/e->a;
  233.     }
  234.     if(y2>pymax) {
  235.         y2 = pymax;
  236.         x2 = (e->c - y2)/e->a;
  237.     }
  238.     if(y2<pymin) {
  239.         y2 = pymin;
  240.         x2 = (e->c - y2)/e->a;
  241.     }
  242.     }
  243.  
  244. /*  fprintf (stderr, "clip_line(): edge %d is (%10.7f, %10.7f) to (%10.7f, %10.7f)\n", numvedges, x1, y1, x2, y2);*/
  245.  
  246.     if (!triangulate)    {
  247.     vedges[numvedges].x1 = x1;
  248.     vedges[numvedges].y1 = y1;
  249.  
  250.     vedges[numvedges].x2 = x2;
  251.     vedges[numvedges].y2 = y2;
  252.  
  253.     vedges[numvedges].nbr1 = (r1 != NULL ? r1->sitenbr : -9998);
  254.     vedges[numvedges].nbr2 = (r2 != NULL ? r2->sitenbr : -9999);
  255.  
  256.     if (r1 != NULL && r2 != NULL)    {
  257.         vedges[numvedges].xm = AVG (r1->coord.x, r2->coord.x);
  258.         vedges[numvedges].ym = AVG (r1->coord.y, r2->coord.y);
  259.         }
  260.  
  261.     if (debug)
  262.         printf ("clip_line puts edge induced by %d and %d\n", r1->sitenbr, r2->sitenbr);
  263.  
  264.     if (numvedges < MAXEDGES)
  265.         numvedges++;
  266.     else    {
  267.         fprintf (stderr, "clip_line(): edge list overflow!\n");
  268.         exit (-1);
  269.         }
  270.     }
  271. }
  272.  
  273. /* load_vsites(): 
  274.     accept the n voronoi sites (x_n, y_n)
  275.     calculate and store the voronoi diagram over the n sites,
  276.     clipping all infinite edges to bbox: [xmin, ymin, xmax, ymax].
  277.     
  278.     note: if (xmin,ymin,xmax,ymax are all == 0.0), OR
  279.     if these do not enclose the data, a bounding box
  280.      will be computed over the input.
  281.  
  282.     returns:
  283.     -1 if error
  284.      0 otherwise
  285. */
  286. int
  287. load_vsites (n, usites, uxmin, uymin, uxmax, uymax)
  288. int n;
  289. float usites[][2];    /* sites in x,y order */
  290. float uxmin, uymin, uxmax, uymax;
  291. {
  292. int k, compute_bbox, sid, tid;
  293. float dx, dy, x, y;
  294.  
  295.     if (n >= MAXSITES)    {
  296.     fprintf (stderr, "load_vsites(): can't handle >= %d sites.\n", MAXSITES);
  297.     return (-1);
  298.     }
  299.  
  300.     compute_bbox = (uxmin == 0.0) && (uymin == 0.0) && (uxmax == 0.0) && (uymax == 0.0);
  301.  
  302.     /* copy the sites into GLsites and set global nsites */
  303.     for (k = 0; k < n; k++) {
  304.     GLsites[k].x = usites[k][0];
  305.     GLsites[k].y = usites[k][1];
  306.     GLsites[k].pid = k;
  307.     }
  308.  
  309.     nsites = n;
  310.  
  311.     /* sort GL sites lexicographically by position */
  312.     sortGLsites();
  313.  
  314.     /* copy sorted GLsites into voronoi alg sites */
  315.     bboxinit();
  316.  
  317.     /* now, if user punted on bbox calculation, OR if user bbox does not truly enclose
  318.        user data, we use our bbox (computed in initbbox).  otherwise we take the user's. */
  319.     if (!(compute_bbox || uxmin > xmin || uymin > ymin || uxmax < xmax || uymax < ymax))    {
  320.     pxmin = uxmin;
  321.     pymin = uymin;
  322.  
  323.     pxmax = uxmax;
  324.     pymax = uymax;
  325.     }
  326.  
  327.     xminymax [0] = xminymin[0] = pxmin;
  328.     xminymin [1] = xmaxymin[1] = pymin;
  329.  
  330.     xmaxymax [0] = xmaxymin[0] = pxmax;
  331.     xmaxymax [1] = xminymax[1] = pymax;
  332.  
  333.     corner[0] = (VERTTHETA *)xminymin;
  334.     corner[1] = (VERTTHETA *)xmaxymin;
  335.     corner[2] = (VERTTHETA *)xmaxymax;
  336.     corner[3] = (VERTTHETA *)xminymax;
  337.  
  338.     /* now: set the floating point tolerance values P*** to be 1 or 2 significant bits inside the p*** */
  339.     /* be careful to use RELATIVE error; that is, to scale the tolerance to the bbox values themselves */
  340.     /* now, if some user puts points way out in left field, our technique handles the ranges correctly */
  341.     {
  342.     float dx = (pxmax - pxmin) * (1.0 / (4096.0));      /* twelve binary digits out */
  343.     float dy = (pymax - pymin) * (1.0 / (4096.0));      /* twelve binary digits out */
  344.  
  345.     Pxmin = pxmin + dx;
  346.     Pxmax = pxmax - dx;
  347.  
  348.     Pymin = pymin + dy;
  349.     Pymax = pymax - dy;
  350.     }
  351.  
  352.     /* compute inverse of external->internal sid mapping */
  353.     for (sid = 0; sid < nsites; sid++)
  354.     inverse[sid] = GLsites[sid].pid;
  355.  
  356.     /* zero list lengths out */
  357.     vdinit ();
  358.  
  359.     /* run the voronoi code, no triangulate */
  360.     triangulate = FALSE;
  361.     voronoi (nextone);
  362.  
  363.     /* RE-copy sorted GLsites into voronoi alg sites */
  364.     bboxinit();
  365.  
  366.     /* run the voronoi code, do triangulate */
  367.     triangulate = TRUE;
  368.     voronoi (nextone);
  369.  
  370.     /* invert the integer id's in sites[], for find_dtriangles() */
  371.     /* and restore the original vertex values (from GLsites)     */
  372.     for (sid = 0; sid < nsites; sid++)    {
  373.         sites[sid].sitenbr = GLsites[sid].pid;
  374.     sites[sid].coord.x = GLsites[sid].x;
  375.     sites[sid].coord.y = GLsites[sid].y;
  376.     }
  377.  
  378.     return (0);
  379. }
  380.  
  381. static VERTTHETA vtlist[1024], slist[1024];
  382. static int  vtnum;
  383.  
  384. int
  385. vtcomp (vt1, vt2)
  386. VERTTHETA *vt1, *vt2;
  387. {
  388.     if (vt1->theta < vt2->theta)
  389.     return (-1);
  390.     else if (vt1->theta > vt2->theta)
  391.     return (1);
  392.     else
  393.     return (0);
  394. }    
  395.  
  396. /*
  397.     find_vregion(sid, plen, pverts)
  398.     given a site id 'sid' from 0..nsites-1 inclusive, 
  399.     returns the voronoi polygon associated with that site
  400.     in the array 'pverts', and its length on call stack.
  401.  
  402.     the vertices are returned in counterclockwise order.
  403.  
  404.     returns:
  405.     -1 if error condition
  406.     plen > 2 [i.e., the # of verts in the voronoi polygon] otherwise
  407. */
  408. int
  409. find_vregion (vsid, pverts)
  410. int vsid;
  411. float pverts[][2];
  412. {
  413. int sid, b, k, vnum, bnd1, bnd2, bdiff, sleft, lag, lead, p1a, p1b, p2a, p2b;
  414. float x1, y1, x2, y2, theta1, theta2, lasttheta, dtheta, h;
  415.  
  416. /* note that PUTV(v,sid,vid) has the side effect of incrementing vid */
  417. #define PUTV(v,sid,vid)    { pverts[vid][0] = (v)->x; pverts[vid][1] = (v)->y; vid++; }
  418.  
  419.     if (vsid < 0 || vsid >= nsites) {
  420.     fprintf (stderr, "find_vregion(%d) called with illegal site id.\n", vsid);
  421.     return (-1);
  422.     }
  423.  
  424.     /* first thing is to convert user's 'virtual' site id to 
  425.        a 'physical' site id, namely, an index into GLsites */
  426.     for (sid = 0; sid < nsites; sid++)
  427.     if (GLsites[sid].pid == vsid)
  428.         break;
  429.  
  430.     if (sid == nsites)    {
  431.     fprintf (stderr, "find_vregion(%d) can't find requested site id.\n", vsid);
  432.     return (-1);
  433.     }
  434.  
  435.     for (k = 0; k < 4; k++)
  436.     corner[k]->theta = atan2 (corner[k]->y - GLsites[sid].y, corner[k]->x - GLsites[sid].x);
  437.  
  438.     for (vtnum = 0, k = 0; k < numvedges; k++)
  439.     if (vedges[k].nbr1 == sid || vedges[k].nbr2 == sid)   {
  440.         /* add both ends of the edge, and their thetas, to unsorted list (parent is edge k) */
  441.         slist[vtnum  ].e1 = slist[vtnum  ].e2 = k;
  442.         slist[vtnum  ].theta = atan2 (vedges[k].y1 - GLsites[sid].y, vedges[k].x1 - GLsites[sid].x);
  443.         slist[vtnum  ].x = vedges[k].x1;
  444.         slist[vtnum++].y = vedges[k].y1;
  445.  
  446.         slist[vtnum  ].e1 = slist[vtnum  ].e2 = k;
  447.             slist[vtnum  ].theta = atan2 (vedges[k].y2 - GLsites[sid].y, vedges[k].x2 - GLsites[sid].x);
  448.         slist[vtnum  ].x = vedges[k].x2;
  449.         slist[vtnum++].y = vedges[k].y2;
  450.         }
  451.  
  452.     /* now we have a list of vtnum thetas and vertices.  sort it on theta */
  453.     qsort ((char *)slist, vtnum, sizeof (VERTTHETA), vtcomp);
  454.  
  455.     /* next, install the unique slist entries into vtlist */
  456.     lag = lead = 0;
  457.     vtlist[lag] = slist[lead];
  458.     lasttheta = -10.0;
  459.  
  460.     while (lead < vtnum)    {
  461.         if (fabs (slist[lead].theta - lasttheta) > 1E-4) {
  462.         lasttheta = slist[lead].theta;
  463.         vtlist[lag++] = slist[lead++];
  464.         }
  465.     else    {
  466.         vtlist[lag-1].e2 = slist[lead++].e1;
  467.         }
  468.     }
  469.  
  470.     vtnum = lag;
  471. /*
  472.     printf ("\n");
  473.     for (k = 0; k < vtnum; k++)
  474.         printf ("vtlist[%d]: x,y %f,%f; theta %f; edges %d %d\n",
  475.         k, vtlist[k].x, vtlist[k].y, vtlist[k].theta, vtlist[k].e1, vtlist[k].e2);
  476. */
  477.  
  478.     for (vnum = 0, k = 0; k < vtnum; k++)    {
  479.         if (fabs(vtlist[(k + vtnum - 1) % vtnum].theta - vtlist[k].theta) < 1E-4)   {
  480.         fprintf (stderr,"find_vregion(%d): vtlist %d, %d identical?\n", sid, (k + vtnum - 1) % vtnum, k);
  481.         return (-1);
  482.         }
  483.  
  484.         x1 = vtlist[(k + vtnum - 1) % vtnum].x;
  485.         y1 = vtlist[(k + vtnum - 1) % vtnum].y;
  486.     p1a = vtlist[(k + vtnum - 1) % vtnum].e1;
  487.     p1b = vtlist[(k + vtnum - 1) % vtnum].e2;
  488.  
  489.     x2 = vtlist[k].x;
  490.     y2 = vtlist[k].y;
  491.     p2a = vtlist[k].e1;
  492.     p2b = vtlist[k].e2;
  493.  
  494.     /* now: if the two vertices come from the same edge, output and continue */
  495.     if (((p1a == p2a) || (p1a == p2b) || (p1b == p2a) || (p1b == p2b))  && (vtnum > 2))    {
  496.         PUTV (&vtlist[k], sid, vnum);
  497.         continue;
  498.         }
  499.  
  500.     /* otherwise must fill in missing corners between x1,y1 and x2,y2 */
  501.     bnd1 = NBOUND(x1,y1);
  502.     bnd2 = NBOUND(x2,y2);
  503.  
  504.     /* find number of CLOCKWISE steps around bbox */
  505.     if (bnd1 >= 0 && bnd2 >= 0)
  506.         bdiff = ((bnd2 - bnd1) + 4) % 4;    /* from 0 to 3 */
  507.     else
  508.         bdiff = 0;
  509.  
  510.     /* if they were on the same bounding box edge, output and continue */
  511.     if (bdiff == 0) {            
  512.         PUTV (&vtlist[k], sid, vnum);
  513.         continue;
  514.         }
  515.  
  516.     /* special case: exactly two vertices */
  517.     if (vtnum == 2) {
  518.         theta1 = vtlist[0].theta;
  519.         theta2 = vtlist[1].theta;
  520.         dtheta = theta2 - theta1;    
  521.  
  522.         /* add all corners OUTSIDE the angular range of the edge as seen from the site */
  523.         for (b = 0; b < 4; b++)    {
  524.         if ((dtheta >= M_PI) && (corner[b]->theta > theta1) && (corner[b]->theta < theta2))
  525.             vtlist[vtnum++] = *corner[b];
  526.         else if ((dtheta < M_PI) && ((corner[b]->theta < theta1) || (corner[b]->theta > theta2)))
  527.             vtlist[vtnum++] = *corner[b];
  528.         }
  529.  
  530.         /* resort the (small) vertex list by theta */
  531.         qsort ((char *)vtlist, vtnum, sizeof (VERTTHETA), vtcomp);
  532.  
  533.         /* and output */
  534.         for (b = 0; b < vtnum; b++)
  535.         PUTV (&vtlist[b], sid, vnum);
  536.  
  537.         break;
  538.         }
  539.  
  540.     for (b = 0; b < bdiff; b++)
  541.         PUTV (corner[(bnd1 + b) % 4], sid, vnum);
  542.  
  543.     PUTV (&vtlist[k], sid, vnum);
  544.     }    /* k 0 ..vtnum */
  545.  
  546.     return (vnum);
  547. }
  548.  
  549. /*
  550.     int
  551.     find_dtriangles (**dtris)
  552.  
  553.     returns:
  554.     -1 if error condition, *dtris == NULL
  555.     o/wise, # of delaunay triangles, *dtris == array of TRIS (see voronoi.h)
  556. */
  557. int
  558. find_dtriangles (dtris)
  559. TRI **dtris;
  560. {
  561.     if (numtris <= 0)    {
  562.     *dtris = NULL;
  563.     return (-1);
  564.     }
  565.     else    {
  566.     *dtris = tris;
  567.     return (numtris);
  568.     }
  569. }
  570.